/*
 *  Matrix.cpp
 *  T3nsors
 *
 *  Created by Michael Barriault on 10-09-23.
 *  Copyright 2010 MikBarr. All rights reserved.
 *
 */

#include <T3nsors/Matrix.h>
#include <T3nsors/Macros.h>
#include <T3nsors/Others.h>
#include <cstdarg>

T3::Matrix::Matrix(void) {
	return;
}

T3::Matrix::Matrix(int iN, std::string iid) {
	*this = T3::Matrix(Domain(1,0,0.,0.),iN,iid);
}

T3::Matrix::Matrix(T3::Domain iO, int iN, std::string iid) {
	O = iO;
	N = iN;
	id = iid;
	T.assign(N, Vector(O,N,id));
	FOR(n,N)
		T[n].id = T[n].id + itos(n);
}

T3::Vector& T3::Matrix::operator()(int i) {
	return T[i];
}

T3::Scalar& T3::Matrix::operator()(int i, int j) {
	return T[i](j);
}

real& T3::Matrix::operator()(int m, int n, int a, ...) {
	va_list va;
	va_start(va, a);
	int t=0;
	FOR(i,O.N) {
		int A = a;
		for ( int j=i-1; j>=0; j-- )
			A *= O(j);
		t += A;
		a = va_arg(va, int);
	}
	va_end(va);
	return T[m](n)[t];
}

T3::Matrix& T3::Matrix::operator+=(T3::Matrix& x) {
	CheckDims(O, x.O); CheckDims(N, x.N);
	PFOR(a,x.N) T.at(a) += x(a);
	return *this;
}

T3::Matrix diff(T3::Vector x) {
	T3::Matrix dx(x.O, x.O.N, (std::string)"d"+x.id);
	PFOR(p,dx.N)
		dx(p) = diff(x(p));
	return dx;
}

T3::Scalar div(T3::Vector x) {
	T3::Matrix Dx = diff(x);
	T3::Scalar dx(x.O, (std::string)"div"+x.id);
	FOR(o,x.O()) FOR(a,x.N) dx[o] += Dx(a,a)[o];
	return dx;
}

T3::Vector curl(T3::Vector x) {
	T3::Matrix Dx = diff(x);
	T3::Vector dx = T3::Vector(x.O, x.N, (std::string)"curl"+x.id);
	dx(0) = Dx(2,1)-Dx(1,2);
	dx(1) = Dx(0,2)-Dx(2,0);
	dx(2) = Dx(1,0)-Dx(0,1);
	return dx;
}

T3::Matrix LinCom(real a, T3::Matrix w, real b, T3::Matrix x, std::string id) {
	T3::Domain O = Merge(w.O,x.O);
	T3::Matrix v(O, w.N, id);
	PFOR(n,v.N) FOR(m,v(n).N) FORO v(n,m)[o] = a*w(n,m)[o]+b*x(n,m)[o];
	return v;
}

T3::Matrix LinCom(real a, T3::Matrix w, real b, T3::Matrix x, real c, T3::Matrix y, std::string id) {
	CheckDims(w.N, x.N); CheckDims(x.N, y.N);
	T3::Domain O = Merge(Merge(w.O,x.O),y.O);
	T3::Matrix v(O, w.N, id);
	PFOR(n,v.N) FOR(m,v(n).N) FORO v(n,m)[o] = a*w(n,m)[o]+b*x(n,m)[o]+c*y(n,m)[o];
	return v;
}

T3::Matrix LinCom(real a, T3::Matrix w, real b, T3::Matrix x, real c, T3::Matrix y, real d, T3::Matrix z, std::string id) {
	CheckDims(w.N, x.N); CheckDims(x.N, y.N); CheckDims(y.N, z.N);
	T3::Domain O = Merge(Merge(Merge(w.O,x.O),y.O),z.O);
	T3::Matrix v(O, w.N, id);
	PFOR(n,v.N) FOR(m,v(n).N) FORO v(n,m)[o] = a*w(n,m)[o]+b*x(n,m)[o]+c*y(n,m)[o]+d*z(n,m)[o];
	return v;
}
